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We investigate the relative phase coherence properties and the occurrence of demixing insta- 
bilities for two mutually interacting and time evolving Bose-Einstein condensates in traps. Our 
treatment naturally includes the additional decoherence effect due to fluctuations in the total num- 
ber of particles. Analytical results are presented for the breathe-together solution, an extension of 
previously known scaling solution to the case of a binary mixture of condensates. When the three 
coupling constants describing the elastic interactions among the atoms in the two states are close 
to each other, a dramatic increase of the phase coherence time is predicted. Numerical results are 
presented for the parameters of the recent JILA experiments. 
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Since the recent experimental observation of Bose-Einstein condensation in dilute atomic gases j^] , much interest 
has been raised about the coherence properties of the condensates. Considerable attention has been devoted to the 
^ matter of the relative phase between two condensates: how this phase manifests itself in an interference experiment 
, how it can be established by measurement |^J^ , and how it evolves in presence of atomic interactions 1^-0] and 
in presence of particle losses ^. 
,; As it was proved in recent experiments performed at JILA, binary mixtures of condensates represent an ideal system 
' to study the phase coherence properties of Bose-Einstein condensates ||^ . In these experiments two condensates in 
0^ , two different internal atomic states are created with a well defined relative phase. After a time r during which the 
' condensates evolve in the trapping potentials, one mixes coherently the two internal atomic states which makes the 
, two condensates interfere; from the spatial interference pattern one gets the relative phase of the two condensates. By 
repeating the whole experimental process, one has access to the distribution of the relative phase after an evolution 
^ , time r, so that one can investigate phase decoherence as function of time. 
'■^ ' The interaction between the two condensates in the JILA experiment gives rise to a rich spatial separation dynamics 
^1 , between the two condensates [|l^, which complicates the theoretical study of the relative phase dynamics. As a 
consequence previous theoretical treatments of the phase decoherence processes, dealing essentially with steady state 
condensates, as in fill , cannot a priori be applied to the experimental situation. 
. ^ ' ^ treatment of the phase coherence of two interacting, non stationary, condensates can be found in |]l2| , with two 
, important differences as compared to the present situation of interest: (1) In jl^ the condensates are subject to a 
■ continuous coherent coupling of amplitude A; results are obtained from a perturbative expansion in powers of 1/A 
5^ , and cannot be simply extended to the present A = case. (2) In fl^ all the coupling constants gaa: gab, gbb between 
the two internal atomic states a and b are assumed to be equal. 

In this paper we propose a formalism to study the relative phase dynamics of interacting and dynamically evolving 
Bose-Einstein condensates initially at zero temperature. 

We present the general method in section^. It consists in expanding the initial state on Fock states, and in evolving 
each Fock state in the Hartree-Fock approximation. We determine the time dependence of the phase collapse for a 
binary mixture of condensates, due to (1) fluctuations in the relative number of particles between the condensates, 
intrinsic to the initial state with well defined relative phase, and (2) fluctuations in the total number of particles. In 
the next two sections we apply this general formalism to two limiting cases that can be treated analytically. 



The first case, in section [II , considers a particular solution of the non- linear Schrodinger equations for the con- 
densates wavef unctions; in this solution the two condensates remain spatially superimposed as they breathe in phase, 
provided that dynamical stability conditions (that we determine) are satisfied. We find that phase decoherence can be 
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highly reduced with respect to non mutually interacting condensates when the three coupling constants gaa, 9ab, gtb 
between atoms in the two internal states a, b are close to each other. 

In section [V, we therefore study in a more general case (not restricted to the breathe-together solution) the 
dynamics of the relative phase for a mixture of condensates for close coupling constants. Our treatment requires also 
in this case the absence of demixing instability, a point that we discuss in detail. 

Finally we discuss the case of the JILA experiment in section ^ This case, that corresponds to close coupling 
constants in a regime of demixing instability, is more difficult to analyze. The predicted phase collapse time depends 
on the fluctuations of the total number of particles ; it is on the order of 0.4 second for Gaussian fluctuations of 8%. 



II. GENERAL METHOD 



In this section, we introduce a gedanken experiment to characterize phase coherence between two condensates: 

the relevant quantity is the interference term {^p^, {r,t)'ipair,t)) between the atomic fields of the two condensates a 
and b. Subsequently we express this interference term in the Hartree-Fock approximation, assuming an initially well 
defined relative phase between the condensates. After a further approximation on the modulus and the phase of the 
condensate wavefunctions, we determine the decay with time of the interference term due to atomic interactions; we 
arrive at the simple results Eq. ( |l8| ) for a fixed total number of particles and Eq. (|2^ ) for Gaussian fluctuations in the 
total number of particles. 



A. Considered gedanken experiment 

The experimental procedure we consider to measure the phase coherence is inspired by recent experiments at JILA 
A condensate is first created in a trap in some internal atomic state a; the corresponding condensate wavefunction 
in the zero temperature mean-field approximation is (pQ, a, stationary solution of the Gross-Pitaevskii equation: 

^2 

Af^o = -^^'^o + [Ua{r) + Ngaa\M^]<l>0- (1) 

In this equation N is the number of particles, gaa is the coupling constant between the atoms in the internal state a, 
related to the scattering length aaa by gaa = 4:711%^ aaa/m; Ua is the trapping potential seen by the atoms in a and /i 
is the chemical potential. Note that we have normalized 0o to unity. 

At time t = a resonant electromagnetic pulse transfers in a coherent way part of the atoms to a second internal 
state b. The state of the system is then given in the Hartree-Fock approximation by 

|vl'(0)> = [ca|a,(/.o)+C6|6,<^o)]'^ (2) 

with |cap -|- |cbp = 1. As we assume a Rabi coupling between a and b much more intense than ^/Ti the atomic 
interactions have a negligible effect during the transfer so that the amplitudes Ca,b depend only on the pulse parameters, 
not on the number N of particles. In the iV-particlc state Eq.(E|) the condensate in state a and the condensate in 



state b have a well defined relative phase; we therefore call this state a phase state, in analogy with |13|. 

The two condensates evolve freely in their trapping potentials during the time r. During this evolution we assume 
that there is no coherent coupling between a and b to lock the relative phase of the condensates; in particular the 
only considered interactions between the particles are elastic, of the type a + a a -\- a (coupling constant gaa > 0), 
a + b ^ a + b (coupling constant gab >0),b + b^b + b (couphng constant g^b > 0). We therefore expect a collapse 
of the relative phase for sufficiently long times, due to atomic interactions. 

To test the phase coherence at time r, a second electromagnetic pulse is applied to mix the internal states a and 
b. We assume that this second pulse is a 7r/2 pulse, so that the atomic field operators in the Heisenberg picture are 
transformed according to 

Mr+)^^Mr-) + ^Mr-) (3) 
Mr+) = -^Mr') + ^Mr-), (4) 
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i5 being an adjustable phase. One then measures the mean spatial density pa in the internal state a, averaging over 
many realizations of the whole experiment: 

Pa = mr+)Ur+)). (5) 

The signature of a phase coherence between the two condensates at time r is the dependence of the mean density 
Pa on the adjustable phase 5. More precisely we define the contrast 

^ _ maxgpa - min^pa ^ 2| (-0|(T~)V'a(T~)) | 
m-AXsPa + minspa E5=a,h('/'^(^")^e(T-)) 

The contrast involves the interference term {'4'l{T^)ipa{T^)) which carries the information about the relative phase 
between the two condensates. 



B. Approximate evolution of an initial phase state 



The time evolution in the phase state representation is not simple, as an initial phase state is mapped onto a 
superposition of phase states. It is more convenient to introduce Fock states, that is states with a well defined number 
of particles in a and in b, these numbers being preserved by the time evolution. 

We therefore expand the initial phase state over the Fock states: 



N 



\m) = E 



iV! 



1/2 



c^cf^|7Va:</>o,^ffc:0o) 



(7) 



where we set Ni, = N — Na ■ 

By calculating the evolution of each Fock state in the simplest Hartree-Fock approximation, we get the following 
mapping: 



\Na : ^o,Nb : (t>o) ^ e-*^(^-^''*)/'*| A^a : M^a, Nb;t), Nb : MNa,Nb;t)) 
where the condensates wavefunctions evolve according to the coupled Gross-Pitaevskii equations: 



ihdt(j)e = 

(where e' ^ e) with the initial conditions 



— A + C/,(f) + 7V,g, 



and where the time dependent phase factor A solves: 



dt 



A{Na,Nb;t) 



1 



1 



K9aa / dr\<j)aV--Nigbb / dr\^br-NaNbgab / dr\<iyaV\M 



(8) 



(9) 



(10) 



(11) 



Equation ( |ll|) is derived in Appendix Physically dA/dt is simply the opposite of the mean interaction energy 
between the particles in the Fock state. In the case where the Fock state is a steady state, the need for the phase 
factor A additional to the Gross-Pitaevskii equation is obvious; the exact phase factor is indeed e"'^*/'', where E is the 
energy of the Fock state, whereas the phase factor obtained from the Gross-Pitaevskii evolution is 6"*^^"^"+^''''''^*/'', 
where pa.b is the chemical potential in a, b. 

Using the evolution of the Fock states, and other approximations valid in the limit of large numbers of particles (as 
detailed in the Appendix |^) we obtain for the interference term between the condensates with a well-defined total 
number N of particles: 

N 

{'4'l'4>a)N = CaCl E 



{Na - l)\Nb\ 



\Ca\'^''-^%b?'''4>a{Na,Nb)4>l{Na - l,iV,+ 1) 



(12) 



where Nb ^ N — Na- The exact computation of this sum remains a formidable task, since it involves in principle the 
solution of N different sets of two coupled Gross-Pitaevskii equations. We introduce some simplifying approximations 
in the next subsection. 
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C. Phase collapse for a mixture 



In the present experiments the total number of particles fluctuates from one realization to the other, so that Eq.(p^) 
has to be averaged over N. We assume that the fluctuations of the total number of particles have a standard deviation 
AN much smaller than the mean total particle number TV. As the distributions of the number of particles in a and 
b in a phase state have also a width much smaller than N (typically on the order of N^^^) we can assume than the 
number of particles in a and in b are very close to their average values — \c^\'^N. We now take advantage of this 
property to simplify Eq.(p^). 

We split the condensate wavefunctions in a modulus and a phase 0^; we assume that the variation of the modulus 
can be neglected over the distribution of Na^b, and that the variation of the phase can be approximated by a linear 
expansion around Nf,. We thus get the approximate form for the condensate wavefunctions: 



(j)eiNa,Nb) ~ exp 



J2 iN,,-N,,)idN^,e,)iNa,Ni,) 



(13) 



where 0^ = <l>e{Na = iV„, TVfc = m). 

To this level of approximation the mean densities in the internal states a, b are simply given by 

{i;li,,)N^N,\<P,\^ (14) 

whereas the interference term between the condensates is: 

{^l4'a)N ^ NCaClcfa^b* CXp {z[(iV - N)Xs ~ iV(|c,|2 - |cfc H^rf] } 6^^" 

[\ca\^e'^^ + \cb\^e-'^'f'\ (15) 
In this last expression we have introduced the time and position dependent quantities 

Xs^l [{dN^ + dN,) i0a - Ob)] (iVa, Nb) (16) 

Xd = l [{dN^ - On,) i0a - 0b)] (iVa, Nb). (17) 



The phase xo — hi^^N^ ~ dN^){0a + 0b){Na, Nb) in Eq.(15) is less important as contrarily to Xs,d it is not multiplied 



by N. At time t = all the x's vanish. In the large N limit, the x's are expected to be on the order of flt/hN. 

The factor responsible for the collapse of the contrast at a fixed value of N is the second line of Eq. (|5|) , the factors 
in the first line being of modulus one. As N is large a small variation of Xd from its initial value Xd{t = 0) = is 
sufficient to destroy the interference term. Over the range of the collapse we can therefore expand the exponential of 
±ixd to second order in Xd, obtaining: 

{^l^fa)N ^ Nca4<Pa'Pb*e^p{t[{N~N)[xs + {\caf ~ \cbf)xd]} cxp [-2Nxl\ca\^\cb\^] ■ (18) 

The second exponential factor in this expression allows to determine the collapse time t^^ for a fixed number of 
particles, through the condition 

m\cM\lit^,n ^ 1 (19) 

such that the modulus of the interference term is reduced by a factor e~^/^ from its initial value. The first exponential 
factor in Eq.(|l^ accounts for the phase difference of the interference term for N particles and N particles, as shown 
by the identity: 

Xs + {\Ca\' - \cb\')Xd = [{Oa ~ 6b){N\Ca\\N\cb\'')]^^f, . (20) 

This phase factor can also be understood as a consequence of a drift of the relative phase between two condensates 
at a velocity v{N) depending on the total number of particles: 

yi^N) = ^^^i^ + {N-N) [xs + (|caP - IcHxd] . (21) 

As we shall see in the next subsection fiuctuations in the total number of particles N result in fluctuations of this 
phase factor, providing an additional source of smearing of the phase, as already emphasized in H. 
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D. Effect of fluctuations in the total number of particles 



The effect on the phase cohapse of fluctuations in the total number of particles is obtained by averaging Eq. ( p^ ) 
over the probability distribution of N. To be specific we assume a Gaussian distribution for N. The average can be 
calculated by replacing the discrete sum over N by an integral; we neglect a term proportional to {ANxd)"^ scaling 
as {AN/N)'^ at the collapse time t^^; the resulting modulus of the interference term reads: 



N\CaCl4)a(t>b*\e'X-P ■ 



exp [-2Nxl\ca 



(22) 



N=N 



The first exponential factor in this expression represents the damping of the interference term due to fluctuations in 
the total number of particles; the second exponential factor, already present in Eq.( |lq), gives the damping due to 
fluctuations in the relative number of particles between a and b, as can be seen in Eq.(|l7|). 



E. The steady state case and comparison with previous treatments 

Our treatment can be easily adapted to the case of two initially different condensate wavefunctions (l)a{t — 0) and 

(f)b{t ^ 0). In the particular case of condensates in stationary states, the formulas for the interference term {'ipb'ipa) 
remain the same, and one has 9^ — —^^i^a, Nb)t/h. We can give in this case the explicit expression for the collapse 



time t^^ defined in Eq.(19), assuming a fixed total number of atoms N = N: 

= h [N^/\aCb\\{dN^ - dNjifia " Mh)]] ■ (23) 

For the particular case of non mutually interacting steady state condensates /i^ depends on N,, only, so that the partial 
derivatives in the denominator of Eq.(^3|) reduce to dfia/dNa + dfib/dNb, and we recover the results of ||,^- 

From Eq.(|23|) we see that what matters physically is the change in the difference between the chemical potentials 
of the two condensates when one transfers one particle from one condensate to the other. For this reason the case 
of mutually interacting condensates with close coupling constants can lead to much larger tcS as compared to the 
case of non-mutually interacting condensates. For example, in the case of the JILA experiment assuming that the 
condensates are in steady state, one finds t^^ ~ 3.1 s; by ignoring the interaction between the condensates (setting by 
hand gab = 0) one obtains the much shorter time ~ 0.25 s. The JILA case is analyzed in more detail in our section 0. 

A similar prediction on the reduction of decoherence due to mutual interactions between the two condensates, in 
trapping potentials with different curvatures, was obtained numerically in 

The treatment in Q considers the absolute phase dynamics of a single condensate (in our formalism Cf, = 0) in a 
coherent state. When the condensate wavefunction is stationary one has simply 9a = —fiat/h. From Eq.(p^) with 
AA^ = iV^/^ (as the coherent state has a Poisson distribution for A^) we then find that the phase of the condensate 
order parameter is damped as exp[— A"((i/ia/<iA^)^i^/2?i^] as in 



III. APPLICATION TO THE BREATHE-TOGETHER SOLUTION 



In this section we consider a particular solution of the coupled Gross-Pitaevskii equations for which an approximate 
scaling solution is available when the chemical potential is much larger than the energy spacing between trap levels, 
the so-called Thomas-Fermi regime. We first give the set of parameters for which this solution, that we call the 
breathe-together solution, exists. We then linearize the Gross-Pitaevskii equations around this solution to determine 
its stability with respect to demixing and to obtain the phase coherence dynamics. 

A. Description of the breathe-together solution 

We now determine the set of parameters such that the coupled Gross-Pitaevskii equations Eq. (^) for 

N, = N, = N\c,\^ (24) 
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have a solution with (f>a{r,t) = (j)b{r,t) = ^{f,t). The general condition is that the effective potential, that is the 
trapping potential plus the mean field potential, seen by the atoms in a and in b should be the same. This condition 
is satisfied when: 



Naga 



Ua{r) - C/b(f) = U{r) 

NbQab = NbQbb + Nagab = ^5- 



The resulting Gross-Pitaevskii equation for the condensate wavefunction common to a and b is then: 



ihdt4> 



- — A + U{r}+Ng\^f 
2m 



(25) 
(26) 



(27) 



with the initial condition 0(r, 0) = (f>o[N = N]{r) = (fo, where 4>o is defined in Eq.([T]). 

By rewriting Eq.(p6[) as Na/Nb — [gbb — gab)/{gaa — gab) we see that this equality can be satisfied by choosing 
properly the mixing angle between a and b provided that 



gab < gaa,gbb or 



gab > gaa, gbb- 



(28) 



As we shall see below, only the first case is relevant here, since the second case corresponds to an unstable solution 
with respect to demixing between a and 6. 



B. Linearization around the breathe-together solution 



The strategy to obtain the quantities Xs,d relevant for the phase dynamics is to calculate in the linear approximation 
the deviations 6(pe between the breathe-together solution 4> and neighboring solutions (p^ for slightly different from 



54>, = MNa + SNa, m + SNb) - MNa, Nb) ■ 
From the definitions Eq.(|l^) and Eq.(17) one indeed realizes that in the limit of small SNa'- 



26Na 



Xd 



2SNa 



(29) 

(30) 
(31) 



where SOa^b are the deviations of the phase of the neighboring solutions from the phase of the breathe-together 
solution: 



9a - set = Im 



(32) 



It turns out that homogeneous rather than inhomogeneous linear equations can be obtained for the deviations 
by introducing the quantities: 



Furthermore a partial decoupling occurs for the linear combinations 

Sips = 6ipa + Sfb 
dipd = Sipa - Sipb- 

The sum Scps obeys a linear equation involving 5ipd as a source term: 



(33) 



(34) 
(35) 



ihdtSifis 



2m 



A + U + 2Ng\^f 



6ips 



+ Ng4>'^Sip; + {Nagaa " Nbgbb){W5ipd + 4>^Sip*a). 



(36) 
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The part of this equation involving 5fs is identical to the one obtained for a single condensate with N particles and 
a coupling constant g. The corresponding modes have minimal frequencies on the order of the trap frequency ld for 
an isotropic harmonic trap [11]. 



The difference Sipd obeys the closed equation: 



2m 



ihdtSifd = 
where we have used the identity: 

Nbigbb 



A + U + Ng\^\'^ 



6(pd 



N 



(gaa + gbb - 2gab){\4'\'^S(fd + <t>^ 5ip*d 



gab) = Naigaa - gab) 



NaNb 

N 



{gaa + gbb - '^gab)- 



(37) 



(38) 



As shown in |15[ minimal eigenfrequencies of this equation can be much smaller than lo\ e.g. when all the coupling 
constants are equal, the minimal eigenfrequencies in a harmonic isotropic trap of frequency lo scale as TuJ^ / fi ^ a; in 
the Thomas-Fermi limit. 

For the derivation of the x's it is sufficient to calculate S^pd- The relative phase between the two condensates for 
the considered neighboring solution with — + SN^ particles in the state e is in fact given by: 



1 

2i 



(39) 



as can be checked from the definition Eq. (p 



C. Approximate equations of evolution in the Thomas- Fermi limit 



In the remaining part of this section we assimre an isotropic harmonic trapping potential U{r) = and we 

restrict to the Thomas-Fermi limit fi 3> hcu. 

In the Thomas-Fermi limit it is known Jl6| , p^ that most of the time dependence of the wavefunction if) can be 
absorbed by a time dependent gauge and scaling transform; here we apply this transform to both and dipd- 



p-iv(t) . 



Sifdir, t) 



A3/2(i) 



A3/2(i) 

The scaling factor \{t) solves the second order Newton- type differential equation 

gaa X* 



V g ^ 2x 

A — — r ~ UJ \ 



(40) 
(41) 

(42) 



with the initial condition A(0) — 1, A(0) = 0. The "force" seen by A in Eq.(^2[) derives from the sum of an expelling l/A"^ 
potential due to repulsive interactions between atoms and an attractive A^)otential due to the harmonic confinement 
of the atoms. It leads to periodic oscillations of A, that is to a periodic breathing of the condensates. We have also 
introduced a phase factor involving the time dependent function rj such that fj = p,g / {gaa/^^K). 

In the Appendix ^ we derive approximate evolution equations for and Sipd] we give here only the result. To 
lowest order in the Thomas-Fermi approximation cf) does not evolve and can be approximated by the Thomas-Fermi 
approximation for ^q: 



(f,t)~0o(f) 



15 



1/2 r 



Rl 



1/2 



(43) 



with a Thomas-Fermi radius i?o = ^ 2[i j muP- . The approximate evolution for bipd is conveniently expressed in terms 
of the real function a. and the purely imaginary function /3: 

~* _ - ~ * 
a = <f) 5ipd + (pSifd 

dipd oifd 



f3 



(44) 
(45) 
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These functions have a clear physical meaning. The first one a corresponds to the deviation Spa/ Na — Spt/Nh written 
in the rescaled frame, Sp^ being the deviation of spatial density in the condensate e from the breathe-together solution. 
Apart from a factor i the second function P is the deviation of the relative phase Eq. d39) written in the rescaled frame: 



{Sea-S9b)ir,t)^~iP{r/X,t). (46) 
The equations of evolution for a, (3 are: 

^n^. ( ; ) = m ( I ) (47) 
where the operator L[t) in the Thomas- Fermi approximation reads: 



-At div 



grad(-) 



L{t)^\ " nxA^-T" —J ). (48) 



The initial conditions for a,P at time t — obtained from Eq.(p3[) and Eq.(|T^) are: 



'SNa_6Nb 

f3{0) = 0. (50) 



D. Solution of the Thomas-Fermi evolution equations: stability against demixing 

The strategy to determine the time evolution of a,P is (1) to expand the vector (q;(0), /3(0)) on eigenmodes of the 
operator L{0), and (2) to calculate the time evolution of each eigenmode. 

1. Expansion on modes of L{0) 

Consider an eigenvector {a, (3) of the operator L{0) with the eigenvalue hfl. For one can express the 

component /3 as function of a: 

and obtain the eigenvalue problem for a: 

n a = [-^^ j S[a] (52) 

where we have introduced the Stringari operator: 

S[a] = ^^ISaa g^-^ ^gg^ 



This operator has been studied in |14|. It is an Hermitian and positive operator, with a spectrum quj^ , q non negative 
integer; q is given by 

q = 2n'^ + 2nl + 3n + I (54) 

as function of the radial quantum number n and the angular momentum I. This allows the determination of the 
eigenfrequencies il: 

_ _ 1/2 
/ NaNb {gaa + 9bb - Sffah) \ 1/2 ,rr\ 
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with g > as we have assumed f2 ^ 0. The case of a vanishing ft corresponds to the zero energy mode ao = 0, /?o = 1 
of the operator L(0), as it can be checked from a direct calculation. 

AU the eigenmodes of 1/(0) have been identified. They do not form a complete family of vectors however. The 
vector (a = l,/3 = 0) cannot be expanded on the eigenmodes of L{0). Its first component a is indeed in the kernel 
of the operator S (as S[a] — 0) whereas none of the aq is in the kernel of S {S[aq] = qoj'^aq is not identically zero) 
except when aq is identically zero (for q = {]). The family of eigenvectors of L(0) completed by the additional vector 
(a = = 0) forms a basis. The additional vector is called an anomalous mode, and we set aanom — 1, /Sanom = 0; 
the action of L(0) on the anomalous mode gives the zero energy mode times the constant factor Naigaa — gab) @- 

The mode functions Uq of the operator S are given in jl^] . It turns out that in the expansion of the initial conditions 
for a, (3 Eq.(M3,&G), only the isotropic eigenmodes of L(0) with q = b and the anomalous mode are involved: 



U(o) 



"9=5 



The isotropic eigenmode of S with q = 5, the so-called breathing mode, reads 

aq=5{r) 



Rl 5 



(56) 



(57) 



By Eq.(|5l|) we have /3g=5 = ctq^^NaNblgaa + 9bb — '2gab)/ NhQq^^. For the coefficients of the modal expansion of 
(a(0),/3(0)), we obtain 



Ca 



— . ^anom- 



(58) 
(59) 



2. Evolution of the modes and stability against demixing 

As a second step we determine the time evolution of the modes of the operator L{0). If we consider an eigenstate 
{aq{0), Pq{0)) of L(0) with the eigenenergy hftq and evolve it according to Eq.(|4^), we find that the evolution reduces 
to multiplication by purely time dependent factors Aq{t), Bq{t): 

aqir,t) = Aq{t)air,0) (60) 
Pqir,t) ^ Bqit)pir,0) (61) 



where the factors satisfy the differential equations: 

n 

A2 



tAq - -^Bq (62) 



^Bq = ^Aq (63) 

with the initial conditions Aq{Q) = Bq{0) = 1. Note that the zero energy eigenmode does not evolve, as ilq = 0. The 
anomalous mode has to be integrated separately, leading to 

aa„om(r,i) = l (64) 
/3a„o„.(r, t) = (65) 

We are now able to address the problem of dynamical stability of the breathe-together solution. Dynamical stability 
requires that any small deviation of the <f>^ 's from the breathe-together solution (p should not grow exponentially with 
time. Here an exponential growth of a may correspond to a demixing of the two condensates a and b. 

A first case of instability occurs when gab > gaa,gbb- In this case the eigenfrequencies VLq are purely imaginary 
and Aq , Bq diverge exponentially with time jl^] . We have checked by a numerical integration of the Gross-Pitaevskii 
equations with spherical symmetry that the spatial distribution then acquires a structure of alternating shells of a 
atoms and b atoms (see Fig.^. 
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We suppose from now on that Qab < gaa,gbb- Instability may still occur in this case due to the periodic time 
dependence of the coefficients in the system Eq.(|6^, as shown i n ]20| |. We have studied in more detail the stability 
of the mode q = 5, which is the one populated initially (see Eq.([S6[)); we have found non-zero instability exponents 
a (C5(t) ~ e"*) in a very limited region of the plane (gab/ gaa, gbb/gaa), with very small exponents (cr < 3 x 10~^a;). 
A direct numerical integration of the Gross-Pitaevskii equations did not show any demixing of a and b even at times 
^ (7^^ [|T]. This suggests that the finite instability exponent is an artifact of the Thomas- Fermi approximation. 

We assume in what follows the dynamical stability of the breathe-together solution. 

FIG. 1. Modulus squared of the condensate wavefunctions |<j!>^ i,|(7Va, A'^b) as function of the distance r to the trap cen- 
ter at a time cut ~ 29.5, from a numerical solution of the coupled Gross-Pitaevskii equations in the case of a dynami- 
cally unstable breathe-together solution. We have taken gtb/gaa ~ 1.2 and Qab/gaa ~ 1.5. We have applied a deviation 
SNa = —SNi, = —0.05Na from the exact breathe-together condition. The chemical potential is p. — 2S.9hui. The curve in solid 
line corresponds to (j>a, the dotted curve corresponds to (pt- 



E. Phase dynamics 

In order to calculate the functions Xd,Xs relevant for the relative phase dynamics, we calculate the evolution of 
the deviation Sipd due to a small change in Na, Nt with respect to iVa, Nb, that is we evolve the initial state Eq.js^) 
according to the results of the previous subsection. 

As we assume dynamical stability of the breathe-together solution, the modes with q = 5 perform only oscillations 
in time | |2^ . The relevant contribution for the phase dynamics therefore comes from the anomalous mode, which from 
Eq.(|65|) has a (3 diverging hnearly with time. Assuming /3(r, t) ^ Canom/3anom(i) and using Eqs.(p8|), ( |4^ ) we obtain: 

(SO. - SebW,t) ^ ^^_R^li9aa+g..-2g.b) (SN^ 5Nb_\ f dt_^ ^gg^ 



5 gaa \Na Nb J Jq X^t') 

We specialize this formula with 6Nb — ±SNa and we get from Eq.(|o|), Eq. (|3l|) : 



l_ f d4J_\ gaa + gbb - '^gab /■* dt' 

2h \dNj^^^ gaa Jo A3(i') ^ > 

Xs - {\cb\^ - \Ca\^)Xd. (68) 

We have introduced the derivative of the chemical potential with respect to the total number of particles 
{{dfj./dN){N — N) ~ 2p,/5N in the Thomas-Fermi limit) in order to recover the characteristic time scale for the 
phase collapse of steady state non mutually interacting condensates. Our formula reveals the interest of close cou- 
pling constants, such that gaa + gbb — 'i-gab ^ gaa- In this case Xd is strongly reduced with respect to non mutually 
interacting condensates; A performs small oscillations around the value A = 1 so that the integral over t' can be 
replaced by t. The more general case of close g's not necessarily satisfying the breathe together condition is analyzed 
in the next section. 

We note that the value of Xs as function of Xd could be expected a priori from Eq. (po|) : when Eq.(^) is satisfied, 
the condensate wavefunctions form a breathe-together solution and have therefore a vanishing relative phase for 
Na = N\ca\'^, Nb = A^jchp, whatever the value of iV is. An important consequence is that there is no extra damping 
of the phase coherence due to the fluctuations of the total number of particles (see Eq.(p2|)). 



IV. CASE OF CLOSE COUPLING CONSTANTS 



We consider in this section the case of close coupling constants which leads to a dramatic reduction of the relative 
phase decoherence with respect to the case of non mutually interacting condensates. 

The strategy is to solve approximately the Gross-Pitaevskii equations Eq.(^) for (paiNa, Nb) and (j)biNa, Nb) and 
apply the formulas Eq.(^ 17) directly. For all equal g^s the initial state is indeed a steady state for the Eq.(^ and 
Xs = Xd = 0. For close g's we linearize the Gross-Pitaevskii equations around the initial value in the hydrodynamic 
point of view. 
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A. Linearization in the classical hydrodynamics approximation 

We first rewrite the Gross-Pitaevskii equations Eq.(^) in terms of tlie hydrodynamic variables: 

p,=NMe{Na.Nb)\^ (69) 
V, = - gT&dee{Na,Nb) (70) 

m 

that is densities and velocity fields of the two condensates. We further assume the Thomas-Fermi limit ^ huj and 
neglect the quantum pressure terms as in ||T^ in the time evolution of the velocity fields: 

dtPe + div{p,ve) = (71) 
dtVe + ^ grad v'^ = -— grad [U{f) + p^g^^ + p^'g^e']- (72) 

2 TO 

At this point we introduce the deviations of the densities and velocity fields from their initial values: 

Pe(i) =Pe(0)+<5pe(0 (73) 
V,{t) ^ vM + 5Ve{t) (74) 



where the initial values are given by: 



pS^^)=NMo?{N) (75) 
v,{t = Q)^Q. (76) 

By expanding Eqs.(|7l|),([72|) to first order in the small quantities 5pe, we obtain: 

dt5pe+AiY[NM^5v,]=Q (77) 

dtSve + — grad [5p^g^^ + Sp^^gee'] = -— grad x 

TO TO 

{N,g,, + N,,g,,, -Ngaa). (78) 
By taking the first time derivative of Eq.([77[) we eliminate the velocity field and we get: 

a^^Pe + Mee'S[5p,,] + a, ^ 0. (79) 

e' 

The source terms of these inhomogeneous equations are: 



TO 



div[|(/.or grad |0ol1(A^e5ee + N,,g,,, - Ngaa)- (80) 



The homogeneous part of Eq.(|79|) involves the 2x2 matrix M: 



Ngaa V Ntgab Ntgbb 



and the Stringari operator defined in Eq.(|5^). In order to solve Eq.([79|) we introduce the eigenvectors e± of the matrix 

M with corresponding eigenvalues g±. Consistently with our previous approximations, we calculate, to leading order 
in the differences between the coupling constants, the eigenvalues: 

g+ - 9aa (82) 

9- - ^2 ^9aa + 9bb - '^gab) (83) 
and the components of {Spa,Spb) on the eigenvectors of M : 

Sp+ ~ 6pa + 6pb (84) 

Sp- - -J^Sp<i - J^^Pb- (85) 
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For those hnear combinations we get the decoupled equations: 



9aa 



(86) 



To study the dynamics of the system we expand p± and the source terms a± on the eigenmodes of the Stringari 
operator. It turns out that the source terms are simply proportional to the breathing mode 0^=5 already introduced 
in Eq.([57|). The solution of Eq.(p6[) with the initial conditions 5p± — dtSp± = is then: 



Sp±{r,t) = iV|0o(O)|M±^[l - cos n±t]ag=5{r) 

9± 



with eigenfrequencies and amplitudes given by: 

1/2 



n± = 

A+ = 
A- - 



5g± 

,9a a 

N^gaa + N^gtb + 2NaNbgab 



1 



NaNb 



N^daa 

Nadaa + NhQab " Nbgbb " A^a5afc 



Nga 



(87) 

(88) 
(89) 
(90) 



We note that when the numbers of atoms Na^b satisfy the breathe-together condition Eq.(26) the amplitude A- 
vanishes as expected, since Sp^ = in this case. 



B. Validity of the linear approximation 

In order for our linearized treatment to be valid the deviations Sp± should remain small as compared to the initial 
densities. A first necessary condition to be satisfied is that the eigenfrequencies il± should be real. This imposes the 
positivity of the matrix M, ensured by the positivity of its determinant: 

gib < gaagtb- (9i) 

This condition is known in the case of homogeneous mixtures of condensates as a stability condition against demixing 
psf . To the leading order in the difference between the coupling constants, the condition Eq.(^l]) is equivalent to 

yaa 

+ gbb - 2gab > 0. 

We note at this point that the amplitude A-/g- in the expression for Sp- is a ratio of two small numbers. When 
this ratio is large the system can evolve far from its initial state even in the stable case <?_ > 0: numerical solutions 
of the Gross-Pitaevskii equations confirm this expectation, showing the formation of a crater at the center of one of 
the condensates. We therefore have to impose a second condition: 

|^±^|«1. (92) 

g± 

Finally the present treatment is based on the classical hydrodynamic approximation; by including the quantum pres- 
sure terms in the hydrodynamic equation for the velocity field one can show that this imposes on the eigenfrequencies 

— « (93) 

M 

(see also Appendix This condition can be violated even in the Thomas- Fermi limit, when the eigenvalue almost 
vanishes. In this case one has to include the quantum pressure terms; the decoupling property of 5p± is unaffected; 
for the evolution of Sp^ similar results as in Eq. (p7|) are obtained; we find e.g. Q- ~ 63tuv'^/8p. 
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C. Phase dynamics 



We assume that all the conditions for the validity of the linearized treatment are satisfied so that we can proceed 
to the analysis of the relative phase dynamics. To this aim we write the equation of evolution for the phases 6^ of the 
condensate wavefunctions (j)^ in the classical hydrodynamic approximation: 

dtOe + ^ ( grad = - [C/ + geePe + 9ee' Pe'V^l. (94) 

The equations for the velocity fields previously given are simply the gradient of Eq. ( |9^ ) . By linearizing Eq. ( |9^ ) around 
the initial state 6'e = we obtain for the relative phase: 

ndtiOa - ~ -\(t)oWNagaa + NbQab ~ NbQbb ~ NaQab) (95) 
+ {gab - gaa)5pa + (gbb - gab)Spb- (96) 

The right hand side of this equation is a sum of terms constant in time and of oscillatory functions of time. The 
function 6a — Ob then has two components: an oscillating component and a component diverging linearly with time 
which will dominate for long times. By using the result Eq. (p7|) and the Thomas- Fermi approximation for |0(O)p 
Eq.(^3|) we can calculate the time diverging component and we obtain to leading order in the g's difference: 

" " '^^ [Nagaa- Nbgbb + {Nb- Na)gab]t/h. (97) 



5Ng. 

We now use Eq. (|l7|) and Eq. (po|) to obtain: 



2 \dNJ gaa 



Xs + {\Ca\-'~\cb\')x<i-~T^{^] 



5gaa \dNy 

{\Ca\'^gaa + |c6pffa6 " \cb\^ gbb " \Ca\'^gab) t/h (99) 

where we introduced the derivative of the chemical potential with respect to the total number of particle [dn/ dN){N — 
N) ~ {2/5)p,/N in the Thomas-Fermi limit. As we already found in the particular case of the breathe-together solution 
the constants Xd and Xs governing the relative phase collapse are highly reduced for close g^s with respect to the case 
of non mutually interacting condensates. 



D. Physical interpretation of the results 

We now show that all the previous results of this section can be interpreted in terms of small oscillations of the 
condensates around the steady state. 

Let us introduce the steady state densities pf for the condensates with Na particles in a and Nb particles in b. 
As we are in the case of quasi complete spatial overlap between the two condensates we can use the Thomas-Fermi 
approximation to determine these densities: 

Pa-U ^ pfgaa + pf gab (100) 

Pb-U = p^'^gab + pfgbb (101) 

where p^ are the chemical potentials in steady state. We rewrite these equations in terms of the deviations Spf of 
the steady-state densities from the initial state densities iVe|0oP and in terms of the deviations Sp,, of the chemical 
potentials from p defined in Eq.(^: 

Spa = {Nagaa + Nbgab - iVffaa)|0Op + ^P^gaa + Vb* 5a6 (102) 

5pb = {Nbgbb + Nagab - Ngaa)\M^ + Spagab + Spf gbb- (103) 



Using the fact that the spatial integral of Sp^ vanishes, we get from integration of Eq.(102 103) over the volume of 
\4>o\'^ the approximate relations: 
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2/i 

S^a = _ , , {Nagaa + Nbgab " Ngaa) (104) 
oNgaa 

SfJ-b = TTT—i^bghb + Nagab ~ Ngaa). (105) 

We can therefore check that the relative phase of the condensates in steady state, given by Of — 91^ = —i{S^a~SfJ.b)t/^, 
evolves as in Eq.(|97|). The phase decoherence properties of the evolving mixture are then essentially the same as in 
steady state. 

Moreover we now show that the average (Sp^) of 5pe over the oscillations at frequencies n± coincide with 6pf. 
First, by averaging Eq.([77[) over time we find that the velocity fields have a vanishing time average p4[ . Second, we 
average Eq. (|78|) over time; we find equations for the spatial gradient of (Sp^), which coincide with the spatial gradient 
of Eq.(|lO|, |l^), so that (Spe) = Spf ||. 



V. DISCUSSION OF THE JILA CASE 

In the JILA experiment the values of the three coupling constants between the atoms are known with good precision; 
they are in the ratio [lO| : 

gaa ■■ gab : fffcb = 1.03 : 1 : 0.97. (106) 

No breathe-together solution exists in this case, as gab lies within gaa and gbb- Experimentally half of the particles 
are in the state a so that |cqP — |cf,p — 1/2, and the mean total number of particles is N — 5 x 10^. Although the 
coupling constants are close, the linearized treatment presented in section |^ does not apply either, because condition 
Eq.(p2|) is violated. It is actually found experimentally that the two condensates evolve far from the initial state, with 
formation of a crater in the a condensate while the b condensate becomes more confined at the center of the trap; 



eventually the condensates separate in some random direction |10|. 

To avoid the crater formation and trigger the spatial separation of the two condensates in a reproducible direction a 
small spatial shift is applied to the trapping potential of one of the two states. The two condensates separate, with a 
relative motion exhibiting strongly damped oscillations ||l^ . The system then reaches a steady state that still exhibits 
phase coherence, up to times on the order of 150 ms after the phase state preparation [^. 

A. Time dependent calculations 

We have already studied in ||l5|| the damping of the relative motion between the condensates, by numerical integration 
of the coupled Gross-Pitaevskii equations Eq. (^) . The agreement with the experimental results of |l^] is qualitatively 
good, although the damping in the theory is weaker and incomplete, small oscillations of the condensate wavefunctions 
remaining undamped even at long times. 

We have applied the formalism of section || by numerically integrating the Gross-Pitaevskii equations for the 
parameters of the JILA experiment. The coefficients Xs,Xd^ now complicated functions of time and space, are 
obtained by evolving wavefunctions with slightly different numbers of atoms in a and b. In order to facilitate the 
comparison with the experiments, in which the a:- integrated atomic density Paiu, z) in the internal state a is measured 
after the 7r/2 pulse applied at time r, we calculated the following contrast: 

max^pa - nimsPa 2] J dx (V;t(r-)^b(r-))G'^"^'^| 

CjlLA(y, z) = — : — = —= - - ^ (107) 

maxspa + Uimspa l^e=a,b J "2; iVe 1 0£ | (t ) 

where the interference term Eq.(^) is averaged over a Gaussian distribution of the total number of particles with a 
standard deviation AA^. A direct comparison with the experiment would require the inclusion of the 22 ms ballistic 
expansion, not included in the present numerical calculations. 

Our numerical result for Cjila at the center of the trap for the species a, y = z = 0, is presented in Fig.|^, for 
Gaussian fluctuations in the total number of particles AN/N = 8% corresponding to the JILA experiment [ p6[ , 
together with the pure Gross-Pitaevskii prediction Cgpe obtained by setting all the x's to 0. The Gross-Pitaevskii 
prediction oscillates around (Cqpe) — 0.63. On the contrary the result of the more complete calculation including 
fluctuations in the relative and total number of particles exhibits a damping of the contrast, that we have fitted by 
convenience with the formula Cjila = Cqb"'''*; we obtain Co — (Cqpe) and 7"^ = 0.42s. 
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Note the oscillatory aspect of the curves in Fig.g. More understanding of the structure of the condensate wave- 
functions given by Eq.(^ is required as this point: as detailed in jl^] is a sum of a smooth part, performing 
oscillations with frequencies expected to be close to eigenfrequencies of the steady state condensates and of a 
noisy quasi-stochastic part. The slow oscillatory structure evident on Cgpe comes from this smooth oscillating part 
of the wavefunctions. 

We have also considered the ideal case of a well defined total number of particles. The numerical prediction for the 
contrast Cjila in this case corresponds to a very long lived phase coherence; after a time of 1 second, the contrast is 
still very close to the pure Gross-Pitaevskii prediction. 

FIG. 2. For the parameters of the JILA experiment (not including the 22 ms ballistic expansion), phase contrasts Cjila 
(lower curve) and Cqpe (upper curve) defined in the text, at j/ = z = 0, as function of time in seconds, for the evolving binary 
mixture, with AN = O.OSiV. 



B. Steady state calculations and effect of particle losses 

As the wavefunctions at long times perform mainly oscillations around the steady state we have also tried a much 



simpler steady state calculation (see subsection [IE). During the collapse time the contrast Cjila is a Gaussian in 
time Eq.(p2[), with an initial value 0.958 and with a half-width tc at the relative height e~^/^. We plot in Fig.|| the 
variation of tc as function of the standard deviation AN. As we find Xs/t ~ —7.7 x 10~^ s~^ and Xd/t = —4.5 x 10~^ 
s~^, one has \xs\ — 1x^1/6, so that relatively high values of AA^ are required to observe a significant effect of the 
fiuctuations of the total number of particles on phase decoherence. For AN = O.OSiV the phase decoherence time is 
tc — 0.32 second, close to the result of the time-dependent calculation of Figj^. Note that for such a high value of 
AN/N the decay of the phase contrast in Eq.(p2[) is essentially due to the first exponential factor accounting for the 
smearing of the phase by fluctuations of the total number of particles, the spreading of the phase for a fixed number 
of particles being very small (iVx^(ic)/2 — 0.005). 

We now briefly consider the issue of losses of particles. An intrinsic source of losses in the JILA experiment are 
the inelastic collisions between a atoms and b atoms, resulting in the simultaneous loss of two particles. We estimate 
the mean number (SN) of lost particles from the rate constant K2 for binary inelastic collisions between the states 
\F = l,m = —1) and |-F = 2,m = 2) and from a numerical calculation of the overlap integral / rf^r |0aP|'^6p- 
For the JILA parameters we find (SN) / N = 0.04 at time tc = 0.32 second. One could then naively expect the effect 
of losses on phase coherence to be comparable with the effect of fluctuations of N . 

To test this naive expectation we use the following simple model, inspired by the two- mode model developed in 
and focusing on the effect of the losses on the drift velocity v{N) of the relative phase of the two condensates given in 
Eq. (pT]) . Imagine that the system has initially N condensate atoms and that k binary inelastic collisions have taken 
place at times ti < . . . < tk between time and time t. The shift of the relative phase during t is then given by: 



e=/ dTv{N{T))=v{N)ti+v{N -2){t2-ti) + ...+v{N -2k){t~tk). (108) 

As we do in ||^ we assume a constant mean number of collisions A per unit of time and we average the phase factor e*® 

multiplying the interference term {ipb 4'a) over the probability distribution of the times ti, . . . ,tk and of the number 
of loss events k, 

Pt{tu...,tk;k)^X^e-^' (109) 

to obtain: 

Ke*^')! =exp{-(fc)[l-sin(2x.)/(2x.)]} (110) 



exp 



for Ix.l « 1 (111) 



where 2(fc) = 2Xt = {SN) is the mean number of lost particles during t. At time t = tc = 0.32 second the corresponding 
modulus of the averaged phase factor is on the order of [1 — 4 x 10^^], very close to one: particle losses have a negligible 
effect on the phase coherence at the considered time tc, even if (SN) and AA^ have the same order of magnitude. 



Actually an inspection of the Xs dependent factor in Eq. (E2h and of Eq. (Ill) reveals that these equations have the 



same structure; replacing in Eq.(E2|) the variance AN^ of the total number of particles by the variance Ak^ of the 
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number of loss events (Afc^ = (fc) as k obeys a Poisson law) one recovers Eq.(lll) up to a numerical factor inside the 
exponential. For equally large values of A7V and (fc) the effect of losses on phase coherence is less important than 
that of fluctuations of because Ak^ — (fc) ^ AA^^. 

We have also investigated another source of losses, the collisions of condensate atoms with the background gas 
of the cell. Assuming a lifetime of the particles in the cell of 250 seconds as in we find as well that this loss 
mechanism has a negligible effect on the phase coherence for a time tc — 0.32 second. 

FIG. 3. For the parameters of the JILA experiment (except the 22 ms ballistic expansion), collapse time tc for Cjila b± 
y — z = as sl function of AN/N for zero temperature steady state condensates in the shifted traps. 



VI. CONCLUSION AND PERSPECTIVES 



We have extended previous treatments of the phase dynamics of Bose-Einstein condensates at zero temperature to 
the case of mutually interacting and dynamically evolving binary mixtures of condensates, for a measurement scheme 
of the phase coherence inspired by the JILA experiment. 

We have first applied this extended formalism to the interesting breathe-together solution of the Gross-Pitaevskii 
equations, in which the two condensates oscillate in phase, remaining always exactly spatially superimposed. The 
analytical results for the phase show that a dramatic increase of the phase coherence time can be obtained for close 
coupling constants gaa, gab, gtb describing the elastic interactions between a atoms and b atoms. 

We have also treated analytically the case of close g's, in the absence of demixing instability. Basically the phase 
collapse is identical to the steady state case for the two mutually interacting condensates. 

Finally, we have investigated numerically the more difficult case of JILA. We find a collapse time of the phase on 
the order of 0.4 second, both by a dynamical and a steady state calculation, in the case of Gaussian fluctuations of 
the total number of particles, corresponding to AN/N = 8%. This result for the collapse time is significantly larger 
than the experimental results (no phase coherence measured after 150 ms). We have also estimated in a simple way 
the effect of coUisional losses on phase coherence in the JILA experiment. 

A possible extension of this work could include the effect of the presence of a thermal component in the experiment. 
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APPENDIX A: PHASE CORRECTION TO THE GROSS-PITAEVSKII PREDICTION 



We consider the evolution of the Fock state \Na : (l)a{0),Ni, : <j)b{0)) (with Na.b particles in the internal state a,b). 
The model Hamiltonian we consider contains the one-body Hamiltonians Tig and elastic interactions terms: 



£—a^b 

^gaalpJ ^jj ifalpa + ^^gbb'^Pb i'b i^btpb + gab^'b '^'ai^b (Al) 

where ip^ is the atomic field operator in the internal state e. 

We use the Hartree-Fock type ansatz for the A^-body state vector: 

I*) = e-'^^'y^\Na : Mt), Nb : Mt))- (A2) 
A variational formulation of the Hamiltonian equation 

tnj^\^)^H\^) (A3) 



16 



leads to the Gross-Pitaevskii equations for cj)^ (t) , given in Eq. , up to the undetermined phase factor A corresponding 
formally to a time dependent Lagrange multiplier ensuring the conservation of the norm of j'l'). To determine this 
phase factor A, we multiply Eq.(A3) on the left by the bra (^Pj; we obtain: 



A + ih{Na : Mt),Ni, : Mt)\j^\Na : Mt),Nb : Mt)) = iW^) 



(A4) 



The scalar products are calculated in second quantized formalism, e.g. we find: 



(A5) 



We finally arrive at Eq.(|l l| 



APPENDIX B: DERIVATION OF THE INTERFERENCE TERM 



When the iV-body state vector is initial a phase state Eq. (g) and if one assumes that the Fock states evolve according 
to Eq.(@) one gets the following expression for the interference term between the two condensates: 



TV 



aCl ^ 



iV! 



J[A(N^~l,Nt + l)~A{Na,Nb)]/n 



|2(^a-l)|^, |2Ar, 



Cbr'MNa,NbW,{Na-l,Nb + l) X 



a{Na -l,Nt, + l)\MNa, N^W [{MNa " 1, A^b + l)\MNa, Nf,))] 



(Bl) 



where Ni, = N — Na- In the large N limit, we expand to first order the effect of shifts of iV^ by unity in the last two 
lines of the previous equation: 



MNa - 1, iVb + 1) ~ MNa - 1, Nb) + dN.MNa ~ 1, Nb) 
MNa, Nb) ~ MNa - 1, Nb) + dNAa{Na - 1, Nb) 

A{Na -l,Nb + l)~ A{Na - l,Nb) + dN,A{Na - l,Nb) 
A{Na,Nb) ~ A{Na - l,Nb) + dN^A{Na - l,Nb). 



(B2) 
(B3) 
(B4) 
(B5) 



We then get: 



N 



N = CaCb 



E 



Nl 



N. 

je{N^-i.Nb) 



- (Na-iy.Nbi' 



|2(Afa-l)|^, |2Af, 



Cbr'MNa,Nb)(bUNa~l,Nb + l) X 



where we have introduced the real quantity: 



eiNa,Nb) = -{On, - dNjA{Na,Nb) + i 



J2Ne{UNa,Nb)\{dN, - dNj\MNa,Nb)) 



(B6) 



(B7) 



We calculate the time derivative of Q{Na, Nb) using the Gross-Pitaevskii equations Eq.(||). After lengthy calculations 
we find 



eiNa,Nb) = o. 



(B8) 



In the gedanken experiment considered in this paper, the initial wavefunctions (f>^(t = 0) depend only on Na + Nb so 
that they have a vanishing derivative c^at^ — dN^i a-nd we take initially ^=0; this leads to = 0. The same conclusion 
holds if the initial wavefunctions are real. 
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APPENDIX C: APPROXIMATE EVOLUTION IN THE THOMAS-FERMI LIMIT 



After the gauge and scale transforms Eq. (^ , ^l| ) the equations of evolution for (j> and 6(pd read: 



ihdt4> — 



2toA2 



ihdtSipd = 



2toA2 



9 



9 



9aa>^^ 



U{r)+Ngaa\4>\^-fi 
U{r^+Ngaa\~4>f - 



Sipd 



1 



+ -j^Nbigtb - 9ab){\(t)\ Sipd + (f> S'-Pd)- 



(CI) 



(C2) 



In the Thomas-Fermi limit the terms involving the Laplacian are small; if we neglect them we get for the time 
derivatives of the a and (3 variables defined in Eq.( ^ , ^5|) : 



ihdta — 

ihdtP = -^Naigaa " gab)a- 



(C3) 
(C4) 



The variable a has actually been defined in a way to obtain zero on the right hand side of Eq.(C3). 

The first equation Eq.( |C3| ) is not an acceptable approximation for the evolution of a, we therefore include in dta 
the contribution of the Laplacian terms: 



ihd, 



'tct 



2mA2 



div 



grad (j> grad i 



2|^pgrad/3 



(C5) 



Furthermore, along the lines of reference [0, one can show that (f> has a negligible time evolution in the Thomas-Fermi 

limit; we can then replace (/> by its initial value <j)o and we recover the first line of Eq.(^). 

The second equation Eq.(C4) is an acceptable approximation for the evolution of /3 if the neglected terms, all 

involving spatial derivatives of a, f3 or c^, are small as compared to the right hand side of Eq.(C4), as they are 
expected to be in the Thomas- Fermi limit. Neglecting these terms, we recover the second line of Eq.(|47j). 

In order to estimate the order of magnitude of the neglected terms in the time derivative of f3, we calculate the 
exact derivative: 



1 

A3 



Na{gaa - gab) 



grad P ■ 



2toA2 



grad 4> grad ( 



2toA2| 



Aa — — grad a 



grad <f> grad • 



1 


A0 A^* 




--a 




1 



(C6) 



We replace (jj by <^o- We consider an eigenmode with frequency Qq; from Eq.(|5^) we estimate Aa/a ^ q/Ro- Assuming 
A on the order of 1 we get the condition 



which we can rewrite as 



1 < 



9=1 



(C7) 



(C8) 
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